PHYS123 Group Homework Project

Wolf Mermelstein, Sara Elanchezhian, Mathew Gummere

December 2, 2023

Part G1

Calculate and plot the recursively defined logistics function $x_{n+1} = R x_n (1 - x_n)$.

We do this by creating a generator to yield logistics equation values, and then use it to find specific ranges of logistics equation values for given growth rates. For our samples we use 0.4 for $x_0$.

import numpy as np
import matplotlib.pyplot as plt
from typing import Iterator

"""Generator that continuously yields logistics equation values at x_0 for some growth rate R."""
def compute_logistics_equation_values(R: float, x_0: float) -> Iterator[float]:
    x_prev = x_0
    while True:
        x_prev = R * x_prev * (1 - x_prev)
        yield x_prev
    

"""Plot the logistic equation values for some growth rate R and initial value x_0."""
def plot_logistic_equation(iteration_counts, x_values, R, x_0):
    fig, ax = plt.subplots(figsize=(10, 6))

    # Create a line plot of the iteration counts vs. the logistic equation values.
    ax.plot(iteration_counts, x_values, marker='o', linestyle='-', color='b')
    ax.set_title(f'Logistic Equation (r={R}, x0={x_0})')
    ax.set_xlabel('Iteration')
    ax.set_ylabel('x')
    ax.grid(True)
    
"""
Compute the coordinates for the logistic equation values for some growth rate R and 
initial value x_0.
"""
def compute_logistics_equation_coordinates(R: float, x_0: float, N: int) -> np.ndarray:
    # Compute values using the generator that we previously defined.
    x_values = np.fromiter(compute_logistics_equation_values(R, x_0), float, count=N)
    # Create an array of iteration counts.
    iteration_counts = np.arange(N)
    return iteration_counts, x_values
    
# Plot the logistic equation values for some growth rate R and initial value x_0.
R_values = (1.5, 2, 2.8, 3.2, 3.6)
for R in R_values:
    x_coords, y_coords = compute_logistics_equation_coordinates(R, 0.4, 100)
    plot_logistic_equation(
        x_coords, 
        y_coords, 
        R=R,
        x_0=0.4
    )
    
plt.show()

Part G2

Finding values of $R_n$ where the logisics equation switches from

• Period 1 to period 2, $R_1$
• Period 2 to period 4, $R_2$
• Period 4 to period 8, $R_3$
• Period 8 to period 16, $R_4$
• Period 16 to period 32, $R_5$

"""
Generate bifurcation plot coordinates.

Arguments:
    equilibrium_value_generator: Generator that yields the equilibrium values given some R/G and x_0.
    value_count: Number of values to have for x axis.
    value_min: Initial value for x axis.
    value_max: Maximal value for x axis.
    iterations_to_start_at: Number of values to skip past for each x value.
    
Returns:
    The x and y coordinates.
"""
def compute_bifurcation_plot_coords(
        equilibrium_value_generator,
        value_count, 
        value_min, 
        value_max, 
        x_0=0.4,
        iterations_to_start_at=900, 
        x_per_value_count=100
    ):
    x_coords = np.zeros(value_count * x_per_value_count)
    y_coords = np.zeros(value_count * x_per_value_count)
    next_val = 0

    # For each R value in an evenly distributed range starting at <R_value_min> up to <R_value_max>
    # with R_value_count values overall.
    for x_type_val in np.linspace(value_min, value_max, value_count):
        # Create a generator to compute x values from
        generator = equilibrium_value_generator(x_type_val, x_0)
        
        # Skip past <iteration_to_start_at> number of computed x values
        for i in range(iterations_to_start_at):
            next(generator)
            
        # Generate <x_per_R_value_count> number of x values and add them to the plot at the current
        # R value (recall x axis is R)
        for i in range(x_per_value_count):
            x_coords[next_val] = x_type_val
            y_coords[next_val] = next(generator)
            next_val += 1

    return x_coords, y_coords

"""
Plot some x coordinates against y coordinates over a clipped region.

Args:
    x_coords: The x coords to plot.
    y_coords: The y coords to plot.
    title: The title of the plot.
    x_min: The minimum x cut off.
    x_max: The maximum x cut off.
    width: The width of the plot.
    height: The height of the plot.
    fontsize: The fontsize of the plot.
    dot_size: The size of the dots in the plot.
"""
def plot_diagram(
    x_coords, 
    y_coords, 
    title="",
    x_min=None, 
    x_max=None, 
    width=50, 
    height=30,
    fontsize=24,
    dot_size=0.05,
    gridline_count=30,
):
    # Create the plot
    fig, axes = plt.subplots(figsize=(width, height))

    # Clip the plot if we have a min or max provided
    if x_min or x_max:
        axes.set_xlim(x_min, x_max)

    # Plot the coordinates
    axes.scatter(x_coords, y_coords, marker='o', linestyle='-', color='b', s=dot_size)

    # Set the title and labels of the plot
    axes.set_title(f'{title} (from x={x_min} to x={x_max})', fontsize=fontsize)
    axes.set_xlabel('R')
    axes.set_ylabel('X_eq')
    
    # Specify the positions of our vertical gridlines by disabling the default x axis gridlines
    # and then manually adding in our own.
    axes.grid(axis='y', linestyle='--')
    gridlines = np.linspace(*axes.get_xlim(), gridline_count)
    axes.set_xticks(gridlines)
    for gridline in gridlines:
        axes.axvline(x=gridline, color='black', linestyle='-', linewidth=.3)

# Compute the bifurcation plot coordinates
x_coords, y_coords = compute_bifurcation_plot_coords(
    compute_logistics_equation_values,
    20000, 
    2.5, 
    4,
)

"""An alias to plot a slice of the bifurcation diagram."""""
def plot_logistics_slice(x_min=None, x_max=None):
    plot_diagram(
        x_coords, 
        y_coords, 
        "Bifurcation Diagram", 
        width=300, 
        height=62,
        dot_size=.01,
        gridline_count=100,
        x_min=x_min,
        x_max=x_max,
    )
    
# Plot the bifurcation diagram and a slice of its tail so that later bifurcations are more
# clearly visible.
plot_logistics_slice()
plot_logistics_slice(3.5, 4)
    
plt.show()

Part G3

Repeat the exercise in HW G-2 but instead of using the logistics equation use the "laser" mapping $x_{n+1} = g x_n (1 - \tanh{(x_n)})$. Find the values of $G_n$ where the logistics equation switches from period 1 to period 2.

To achieve this we redefine a new generator, compute_laser_equation_values, to compute values, which we feed into the same x-y coordinate computing function, compute_bifurcation_plot_coords, that we defined earlier. We then plot various subsections to show period switches more clearly.

from numpy import tanh

"""Generator that continuously yields laser equation values at x_0 for some growth rate R."""
def compute_laser_equation_values(G: float, x_0: float) -> Iterator[float]:
    x_prev = x_0
    while True:
        x_prev = G * x_prev * (1 - tanh(x_prev))
        yield x_prev

# Compute the bifurcation plot coordinates for the laser equation
x_coords, y_coords = compute_bifurcation_plot_coords(
    compute_laser_equation_values,
    50000, 
    1, 
    10,
)

# Plot the entire bifurcation diagram
plot_diagram(
    x_coords, 
    y_coords, 
    "Laser Bifurcation Diagram, all periods", 
    width=100, 
    height=40,
    x_min=1,
    x_max=10,
    dot_size=.015
)

"""Plot a slice of the bifurcation diagram for the laser equation."""
def plot_laser_slice(x_min, x_max, title):
    plot_diagram(
        x_coords, 
        y_coords, 
        f"{title} Laser Bifurcation Diagram", 
        width=48, 
        height=26,
        x_min=x_min,
        x_max=x_max,
        dot_size=.04
    )
    
# Plot slices of the bifurcation diagram for the laser equation for different period switches
plot_laser_slice(5, 5.2, "Laser Bifurcation Plot, Period 0-1")
plot_laser_slice(8, 8.2, "Laser Bifurcation Plot, Period 1-2")
plot_laser_slice(9, 9.15, "Laser Bifurcation Plot, Period 2-3")
plot_laser_slice(9.2, 9.38, "Laser Bifurcation Plot, Period 3-4-5")

Part G5

Consider the logisics equation $x_{n+1} = R x_n (1 - x_n)$ with the growth parameter set to $R = 3.6$.

Define a sequence $x_{true}$ to be set of the $x_n$ calculated with $R = 3.6$ and $x_0 = 0.5$

Define a sequence $x_{pred}$ to be set of the $x_n$ calculated with $R = 3.6$ and $x_0 = 0.5 + \epsilon $.

The idea is that even a small error $\epsilon$ in the initial conditon will lead to predicted values that vary wildly from the "true" values.

"""Add an error to an x_0 value."""
def with_error(x_0, error):
    return x_0 + error

"""Build an array of x values for a given error value."""
def compute_x_values(x_0, error=0):
    return np.fromiter(compute_logistics_equation_values(R, with_error(x_0, error)), float, count=N)

"""Superpose an additional plot onto the main plot axes."""
def plot_y_axis(axes, x_values, y_values, color, label, width=.5):
    axes.scatter(x_values, y_values, color=color, label=label, linewidth=width, s=.1)
    axes.set_ylabel(label, color=color)

# Create the plot
fig, axes = plt.subplots(figsize=(40, 16))
plt.title("Effect of x_0 errors on bifurcations")

# Create y axis for the various errors
errors = (0, 0.05, 0.005, 0.0005)
colors = ('black', 'blue', 'red', 'green')

# Plot the data
for error, color in zip(errors, colors):
    x_coords, y_coords = x_coords, y_coords = compute_bifurcation_plot_coords(
        compute_logistics_equation_values,
        3000, 
        2.5, 
        4,
        x_0=with_error(x_0, error),
    )
    label = f'Error = {error}'
    plot_y_axis(axes, x_coords, y_coords, color=color, label=label)
    
# Add a legend
axes.legend(loc='upper left', fontsize=28)
axes.grid(True)

# Show the plot
plt.show()